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When dealing with system-reservoir interactions in an open quantum system, such as a photo¬ 
synthetic light-harvesting complex, approximations are usually made to obtain the dynamics of the 
system. One question immediately arises: how good are these approximations, and in what ways 
can we evaluate them? Here, we propose to use entanglement and a measure of non-Markovianity as 
benchmarks for the deviation of approximate methods from exact results. We apply two frequently- 
used perturbative but non-Markovian approximations to a photosynthetic dimer model and com¬ 
pare their results with that of the numerically-exact hierarchy equation of motion (HEOM). This 
enables us to explore both entanglement and non-Markovianity measures as means to reveal how 
the approximations either overestimate or underestimate memory effects and quantum coherence. 

In addition, we show that both the approximate and exact results suggest that non-Markonivity 
can, counter-intuitively, increase with temperature, and with the coupling to the environment. 


I. INTRODUCTION 

Modelling and understanding the non-equilibrium dy¬ 
namics of open quantum systems is a ubiquitous problem 
in physics, chemistry and biology [1, 2]. In such sys¬ 
tems, the environment is usually composed of huge num¬ 
bers of microscopic constituents, an exact description of 
which is challenging. One can invoke intensive computa¬ 
tional techniques, such as path-integral formalisms [3-6], 
Monte Carlo algorithms [7], the hierarchy equations of 
motion (HEOM) [8-10], the reaction-coordinate method 
[ 1, 12] and others, to explicitly and exactly propagate 
the quantum state of a complete system-environment 
model. However, a common drawback of these exact 
numerical solutions is their demanding computational re¬ 
source requirements, which can scale badly depending on 
the spectral density of the environment being modelled, 
the number of independent baths the system is coupled 
to, or the complexity of the system itself. 

To simplify the problem and gain useful physical in¬ 
sight, approximations are usually made to reduce the 
system dynamics to that of a relatively few degrees of 
freedom. In that regard, much effort has been devoted 
to develop quantum master equations (QMEs) which de¬ 
scribe these reduced degrees of freedom in various limits. 
Redfielcl theory [13] provides one with QME based on (to¬ 
gether with a secular approximation) a second-order per¬ 
turbation approximation in the system-environment cou¬ 
pling. For the strong-coupling limit, Forster theory [14— 
16], adopts a diffusion-rate equation [ 7] to describe the 
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incoherent transport phenomenon. Nevertheless, these 
conventional Markovian QME treatments cannot capture 
the memory effects of the bath. 

In order to take into account the memory effects, many 
attempts at improving these Markovian QMEs have been 
made. The second-order time-convolution (TC2) [ 8] 
equation is known as a chronological-ordering prescrip¬ 
tion [1! ] or time-nonlocal equation [20, 21]. It is a direct 
generalization of Redfield QME without performing the 
Markov and secular approximations. The second-order 
time-local (TL2) equation is another frequently used 
QME, sometimes called a partial-time-ordering prescrip¬ 
tion [19] or time-convolutionless equation. Some works 
suggest that TL2 shows better performance than TC2 at 
numerically approximating exact results [22]. Neverthe¬ 
less, their respective domains of applicability have not 
been thoroughly investigated yet. 

In each QME model (TC2, TL2), certain approxima¬ 
tions and simplifications are introduced to obtain solv¬ 
able equations. To investigate the deviation of each ap¬ 
proximate QME model from the exact results, we first 
compare the explicit dynamics of these two approxima¬ 
tive QMEs with that of the HEOM (which is considered 
to be numerically exact). We focus on the intermediate 
system-environment coupling regime, which has proven 
to be the most challenging and relevant to the dynam¬ 
ics in realistic systems, and the one at which the regime 
of validity of most approximations breaks down. Both 
approximate methods are perturbative in the system- 
bath coupling, but can in principle harbor memory ef¬ 
fects of the environment. To investigate how well they 
capture these effects we then utilize the concept of the 
Clroi-Jamiolkowski isomorphism [23, 24] to encode com¬ 
plete information on the dynamics of the system into 
the entanglement with an ancilla. By comparing the 



2 


time evolution of the entanglement between system and 
ancilla, and an associated measure of non-Markovianity 
[25], one can find out to what extent the memory effects 
and coherence predicted by each approximative QME de¬ 
viates from being numerically exact. Our results suggest 
that entanglement and non-Markovianity provide a use¬ 
ful benchmark for the performance of such approximative 
treatments. 

In performing this analysis we also discuss several in¬ 
teresting physical trends, including a counter-intuitive in¬ 
crease of non-Markovianity with both temperature and 
with the coupling strength to the environment. We 
attribute this increase to an enhancement of system- 
environment correlations when both the coupling and 
temperature are increased. Additionally, evidence from 
other studies [26-28] suggests that non-Markovian envi¬ 
ronments are capable of sustaining quantum coherence. 
The interplay of these factors finally results in the in¬ 
crease of non-Markovianity with both temperature and 
coupling strength that we see in our results. 


II. THE SPIN-BOSON MODEL 
A. The model 

The spin-boson model [3] is one of the most extensively 
studied models of open quantum systems, and is the one 
we employ here. It describes a spinor-like two-state sys¬ 
tem interacting with a bosonic environment. First, let us 
consider this standard model, which can be divided into 
three components 


#tot — -^sys + H e nv + #i n t- (1) 


The system Hamiltonian, Hgys, is written as 

TT _ ^0 - , _ 

tlsys — ~2 ~ Gz + 


( 2 ) 


where Ja x is the coherent-coupling term, which enables 
the tunneling between the two system quantum states, 
labeled as |1) and | —1), with the energy level spacing Hcoq. 
Usually, one adopts the delocalized basis |x+) and |\_) 
(exciton), which is defined by the following eigenvalue 
problem 


+ Jvx') I X±) = ±7^ I X±)i (3) 

with Q = y /uiq + 4J 2 /h 2 . 

The environment, F/env, is usually modelled as a large 
collection of harmonic oscillators 

He nv = ^ (4) 

k 

where (at) is the creation (annihilation) operator of 
the environment mode k with angular frequency Wk- For 


simplicity, a linear system-environment coupling, H- in ^. 
is adopted throughout this work: 

#int ® (^.9kdk + %kdk) , (5) 

k 


where g k is the coupling constant between the environ¬ 
ment mode k and the system. In most physical prob¬ 
lems, the details of the microscopic description of g k 
are not clear, and one usually employs a spectral den¬ 
sity function, J{ui) = ]T] k |.9k| 2 <5(w — Me), to character¬ 
ize the coupling strength via the reorganization energy 
A = f Q J(io)/ojduj. The physical meaning of the spec¬ 
tral density function can be understood as the density 
of states of the environment, weighted by the coupling 
strengths. Moreover the way in which the environment 
modulates the dynamics of the system is described by 
the correlation function 



coth 


/ hu> \ 
\2k B Tj 


cos cut — i sin uit 


du>. 


( 6 ) 

The real part is related to the dissipation process, while 
the imaginary part corresponds to the response function. 

The statistical properties of the entire system can be 
described by the total density matrix p^oD which con¬ 
tains all the degrees of freedom of the system and en¬ 
vironment. If the correlation between the system and 
environment is negligible, the Born approximation can 
be used and the total density matrix can be factorized 
into 


Ptot(^) — Psys(t) <8> Penvj (7) 


where psys(t) describes the dynamics of the system 
and penv = exp —Henv/ksT 1% is the environment 
density matrix in thermal equilibrium at temperature 
T. Here, k B is the Boltzmann constant and Z = 


Tr exp 


—Hemr/k B T 


is the partition function. 


One notes that when wo, J , and A are comparable, 
this makes the conventional perturbative treatment unre¬ 
liable. In the following, we will adopt the two frequently- 
used perturbative but non-Markovian QME formalisms 
discussed in the introduction and compare their results 
with the exact one in the intermediate-coupling regime, 
as they both begin to break down, and investigate ways 
in which to evaluate their accuracy. 


B. Second-order time-convolution equation (TC2) 


For the Hamiltonian defined above, the time evolution 
of the system density matrix pt 0 t(f) under the TC2 ap¬ 
proximation is expressed as 

~Q t P sysCO = -^ Tr env 


h? 


/0 L 


-®int W’ ^int( T )> Psys{ T ) ® Penv dr. (8) 
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The tilde symbol above an operator denotes the interac¬ 
tion picture with respect to Hgys + Heuv- The interac¬ 
tion Hamiltonian in terms of the delocalized basis can be 
expressed as 

\»){v\ 

® (hg\,a\.e lulkt + hg^d k e~ lulkt ^ , (9) 

where A = (p\a z \v), and p, v = X-hX-- Substituting 
Eq. (9) into (8) with the explicit expansion leads to a 
set of simultaneous integrodifferential equations of the 
density matrix elements Pn, v (t) 

— [e~ lUa ’ f,t p aj p(t)] = —iWa,pe- XU ’ a '* it Pa,p{t) 

+ Y1 f fu-A 1 - T ) er ' lu>y " uT P^A T ) dT - ( 10 ) 

One notes that the memory effects are taken into account 
in terms of the convolution of the memory kernel /^(f— 
r). A detailed expression for this kernel is given in the 
Appendix. 

To solve the simultaneous integrodifferential compo¬ 
nents of Eq. (10), we invoke the Laplace transformation 
£{/} := f 0 f(t)e~ st dt , and transform them into a set of 
algebraic equations. After carefully analyzing the prop¬ 
erties of the poles, the conventional residual theorem en¬ 
ables one to accomplish the inverse Laplace transforma¬ 
tion and move back from Laplace space into the time- 
domain. 


C. Second-order time-local equation (TL2) 


In the TL2 formalism, the system is considered to be 
sluggish, hence the bath feedback on the system dynam¬ 
ics can be neglected by approximating psys{t — r) ss 


exp 


iH S y S T/h psys(t) exp - iH S y S T/h 


This assump¬ 


tion is reasonable because it is impossible for a system to 
change its configuration instantaneously. Consequently 
the system density matrix should be pulled out from the 
integral to obtain the following QME 


foPsysit) - -^Tr< 


H 2 


env 


/o L 


-^int(*)> H int (r),p S ys(i) ® penv dr. (11) 


Similarly, substituting Eq. (9) into (11) with the explicit 
expansion leads to a set of simultaneous differential equa¬ 
tions of the density matrix elements p^^it) 


d_ 

At 


[e-^PaAA = - 


VjJ, 


t,pe Pa,pit) 




Pn,v{d)- 


( 12 ) 


The detailed expression of the memory kernel — r) 

is given in the Appendix. It should be emphasized that 
although p Mi „(t) is pulled out from the integral, Eq. (12) 
is capable of predicting a non-Markovian dynamics be¬ 
cause the time integral of h^^it) results in time-varying 
coefficients in front of p^^it). Whether or not such dif¬ 
ferential equations behave non-Markovianly crucially de¬ 
pends on these time-varying coefficients. 


III. COMPARISONS WITH EXACT RESULTS 

To illustrate the differences of the approximations ex¬ 
plicitly, we apply these two QMEs to a photosynthetic 
dimer model, which has attracted considerable interest 
recently [1, 2, 29-34], We employ the Drude-Lorentz 
spectral density function (the over-damped Brownian os¬ 
cillator model) [10, 35], J(uj) = ( 2 A 7 / 7 r) [w/ (w 2 +7 2 )], 
which has been shown to fit the experimental data well 
[29, 30] and has been successfully used for a range of 
theoretical studies of this type of system [31-34], As 
mentioned in the previous section, the reorganization en¬ 
ergy, A, characterizes the coupling strength to the envi¬ 
ronment, while the quantity 7 determines the width of 
the spectral density. These two parameters have consid¬ 
erable influence on the dynamics of the system. 

In Fig. 1, we show the system dynamics given by (a) 
TC2, (b) HEOM, and (c) TL2 with varying A and tem¬ 
perature T. The other parameters are fixed at w 0 = 70 
cm” 1 , J = 100 cm” 1 , and 7 = 50 cm” 1 (7” 1 = 106 fs). 
These parameters are typical in photosynthetic systems. 
The solid curves in each panel denote the populations 
of the |x+) state with temperatures T = 300 K (black), 
250 K (red), and 200 K (blue), respectively. It can be 
seen that, at higher temperatures, the population of the 
|x+) state transfers to the |x_) state faster than at lower 
temperatures, but there is always a crossing so that the 
thermal equilibrium population of the |x+) state is larger 
at higher temperatures. 

For small values of A, the results of the two QME mod¬ 
els show excellent agreement with that of the HEOM, 
indicating that both TC2 and TL2 perform well in the 
weak system-environment coupling regime and that the 
bath memory effect is insignificant at small A. Moreover, 
the result of TC2 completely coincides with that of the 
HEOM for very small couplings. We show the compari¬ 
son between TC2 (solid curve) and HEOM (dot-dashed 
curve) methods in the inset of Fig. 1(a) for A = 5 cm” 1 , 
and T = 250 K. When A is increased, the TC2 population 
results exhibit vigorous beating and produce oscillatory 
curves up to 800 fs, which is absent in the HEOM result. 
We attribute these oscillations to the over-estimation 
of the coherence by TC2. Apart from these beatings, 
the overall magnitude of the population of HEOM is 
quantitatively better approximated by TC2 than TL2. 
The TL2 model yields monotonically-decaying popula¬ 
tion dynamics that tend to reach thermal equilibrium 
too rapidly. This leads to a significant over-estimation of 
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(a) 



(b) 



(c) 



FIG. 1. (Color online) Time evolution of the populations p++(t) (solid) and coherence |p _|_(t)| (dashed) predicted by (a) TC2, 

(b) HEOM, and (c) TL2 for the spin-boson model with different values of A at temperatures T = 300 K (black), 250 K (red), 
and 200 K (blue). The other parameters are u>o = 70 cm -1 , J = 100 cm -1 , and 7 = 50 cm -1 ( 7 -1 = 106 fs). For small A, both 
QMEs yield excellent results, as expected. The inset in (a) shows the results given by TC2 (solid curve) and HEOM (dot-dashed 
curve) for A = 5 cm -1 , and T = 250 K, illustrating how they almost overlap. However, due to over-estimation of the coherence, 
the result calculated from the TC2 method shows a slightly higher beating behavior in the population dynamics. In contrast, 
for large A the population dynamics predicted by the TC2 method is in better agreement with those of the HEOM, whereas 
the populations given by the TL2 method are somewhat sluggish and tend to approach thermal equilibrium a bit faster. 


|p-+(t)| |p-+(t)| |p-+(t)| 
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the population relaxation rate by TL2, especially at large 
A. This over-estimation of the population relaxation rate 
in Redfield theory has been reported previously [36] , and 
here we gain further insight into its origin by comparing 
to the TC2 results. 

The dashed curves in Fig. 1 denote the absolute value 
of the off-diagonal elements of the system density matrix, 
i.e., the coherence between the |x+) and |y_) states. The 
results from the TC2 method manifestly show the over¬ 
estimation of the coherence even if A is small. When A is 
increased, the over-estimation of the coherence becomes 
quite pronounced. On the other hand, the coherence in 
the TL2 model decays more rapidly, leading to the slug¬ 
gish dynamics discussed above. In summary, the coher¬ 
ence dynamics is better approximated by TL2, and the 
TC2 model may fail in approximating the true coherence 
for large A. However, the overall population decay rate 
predicted by the TC2 is generally more correct than that 
of TL2. It is interesting to note that the TL2 model 
yields an exact QME for a pure dephasing spin-boson 
model (i.e. J = 0) [22] while the TC2 model underesti¬ 
mates the pure dephasing rate, which is in line with our 
findings here. 


IV. BENCHMARK OF APPROXIMATIVE 
QMES 

A. Entanglement and non-Markovianity 

In the previous section, we analyzed how the coher¬ 
ence terms of the two approximations are qualitatively 
different from the HEOM exact results. However, those 
comparisons fail in providing an overall intuitive picture 
about which model performs better as they are basis- 
dependent. In other words, it is possible that one model 
may perform better or worse than another depending on 
the bases used. In this section, we apply a measure of 
the non-Markovianity to develop a bases-free benchmark 
which can quantitatively describe the performance of the 
approximate methods. 

Let us consider an isolated ancilla possessing the same 
degrees of freedom of the system and with which the 
system forms a maximally entangled initial state |'F) = 
Ej= x+ ,x- 75 li) ® | j) (see Fig. 2). If the system evolves 
according to a process £t,Oi then the Choi-Jamiolkowski 
isomorphism [23, 24] guarantees that the extended den¬ 
sity matrix 


Psys,anc(i) — (ft,o ®^anc) (|^ r )(^ f |) 

= E^.o(li)(fc|)®li)(fc|- (13) 

j,k 


contains all the necessary information on the dynamics 
of the system, where lane is the identity process acting 
on the ancilla. The entanglement, E(psys,anc), between 


System 



FIG. 2. (Color online) Schematic illustration of the entangle¬ 
ment measure. We consider a system and a copy of it acting as 
a well-isolated ancilla possessing the same degrees of freedom 
of the system. Initially, they form a maximally-entangled 
state psys,anc(0) = |4/)(4/|. Then the system starts to feel 
contact with its environment (denoted by the gray shadow) 
and evolves according to £t,o, whereas the ancilla is kept iso¬ 
lated. 


the system and the ancilla is a physical quantity which 
is typically very sensitive to environmental effects. 

Another related quantity is the degree of non- 
Markovianity, MM. Recently, many efforts have been 
devoted to construct a proper measure of the non- 
Markovianity [25, 37]. Rivas et al. [25] combine the con¬ 
cept of the divisibility of a quantum process [38, 39] and 
the fact that no local completely positive (CP) operation 
[24] can increase the entanglement E between a system 
and its corresponding ancilla 

E[psys,anc] > E [(fsys ® Zanc)(psys,anc )] • (14) 

Consequently, Rivas et al. [25] proposed that the degree 
of non-Markovianity within a given time interval [0, t] can 
be estimated by 


MM = 
where 


^E[(£ T ,o<8Xanc)(|*><*|)] 


dr — AE# 


(15) 


AE t = E p><*|] - E [(St, 0 ® lane) P)(*|)] ■ (16) 

The non-Markovianity of open-system quantum dynam¬ 
ics can be evaluated at many different theoretical levels 
[25, 37, 40-42], and the quantity ATM is an extremely 
strict indicator of non-Markovianity that measures the 
information exchange in time between the system and its 
environment. For MM to have a non-zero value, explicit 
environmental memory effects must be present. 

Here we compare the time evolution of the entan¬ 
glement, Ej, and the corresponding degree of non- 
Markovianity, MM, for the two approximate system- 
bath models and show how they can provide an inte¬ 
grated picture as to what extent their dynamics deviate 
from the exact results. 
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(a) 



(b) 



(c) 



FIG. 3. (Color online) Time evolution of the concurrence calculated by (a) TC2, (b) HEOM, and (c) TL2 for different values 
of A at temperatures T = 300 K (black), 250 K (red), and 200 K (blue). The other parameters are ui o = 70 cm -1 , J = 100 
cm -1 , and 7 = 50 cm ' 1 ( 7 1 = 106 fs). In general, the concurrence will die out faster for larger A and higher temperatures. 
The coherence over-estimation of the TC2 method is manifested by a concurrence revival around 100 fs for larger values of A, 
whereas HEOM and TL2 produce a monotonically-decreasing concurrence. 


B. Evaluating non-Markovianity 

To analyse the behavior of the non-Markovianity in 
each method, in this section we will show how the con¬ 
currence, a well-known measure for bipartite entangle¬ 
ment [43], between system and ancilla evolves in time 
and how the corresponding non-Markovianity [Eq. (15)] 
depends on the physical parameters of the original spin- 
boson model. 

As an explicit visualization of the integrand in Eq. (15), 
in Fig. 3, we apply the measure to (a) TC2, (b) HEOM, 
and (c) TL2 and show the time evolution of the concur¬ 
rence for different values of A at temperatures T = 300 K 
(black), 250 K (red), and 200 K (blue), respectively. The 
other parameters are wq = 70 cm -1 , J = 100 cm -1 , 


and 7 = 50 cm -1 ( 7 -1 = 106 fs). It can be seen 
that, when increasing the temperature and A, the deco¬ 
herence becomes more pronounced. Hence, the concur¬ 
rence will die out earlier for larger A and higher tem¬ 
perature. As shown in Fig. 3(a), except for A = 5 
cm' 1 , which produces monotonically-decreasing concur¬ 
rence, the TC2 model produces oscillatory curves, in 
which a concurrence revival is exhibited around 100 fs 
and results in a finite degree of non-Markovianity (shown 
later). While in Fig. 3(b) and (c), HEOM and TL2 pro¬ 
duce monotonically-decreasing concurrence and generate 
no visible non-Markovianity with this measure. 

In Fig. 4, we show the corresponding measure of the 
non-Markovianity, A fJA, calculated using the time evo¬ 
lution of the concurrence shown in Fig. 3(a). Only TC2, 
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FIG. 4. Non-Markovianity, JVM, obtained by the TC2 
method as a function of temperature. The parameters are 
the same as used in Fig. 3. Among the three methods inves¬ 
tigated in this work, only TC2 at higher A generates non-zero 
non-Markovianity for these parameters. At A = 5 cm -1 TC2 
correctly produces the expected Markovian dynamics for this 
regime. 


for larger A values, leads to non-zero non-Markovianity, 
while TC2 at A = 5 cm -1 , HEOM, and TL2 generate null 
results due to the monotonically-decreasing concurrence. 
This comparison not only shows that the TL2 yields a 
better approximation to the HEOM dynamics, but also 
explicitly demonstrates the degree to which TC2 devi¬ 
ates from HEOM. We again attribute this deviation to 
the over-estimation of coherence shown in Fig. 1. In ad¬ 
dition, it can be seen in Fig. 4 that JVM tends to increase 
with increasing A and temperature. We will investigate 
this below in a regime where the HEOM results exhibit 
similar behavior. 


C. Increase of non-Markovianity with A and 
temperature 

The other two important parameters in our spin-boson 
model are the level spacing uiq and the bath relaxation 
time 7 . The former affects to what extent the state |x+) 
is delocalized, while the latter is related to the correlation 
time of the environment and is directly connected to the 
non-Markovianity of the system. 

In Fig. 5(a), we reduce ojq to 40 cm -1 and fix the other 
parameters at A = 5 cm -1 , 7 = 50 cm -1 , and T = 200 
K. The reduction of uiq leads to a manifest concurrence 
revival around 100 fs in the TC2 concurrence dynam¬ 
ics, a result of stronger delocalization and significant en¬ 
hancement of the coherence effect. In the mean time, the 
concurrence of the HEOM result is still monotonically 
decreasing. The TC2 model further over-estimates this 
enhancement and ends up with finite non-Markovianity 
within all range of temperatures shown in Fig. 5(b). The 
TL2 model predicts almost-Markovian results, besides 
the very small non-Markovianity at low temperatures, 




Temperature (K) 

FIG. 5. (a) The concurrence obtained from the TC2 method 
(dot-dashed), HEOM (dashed), and TL2 (solid), for u>o = 40 
cm” 1 . The reduction of ujo leads to a manifest concurrence 
revival around 100 fs. The concurrence obtained from HEOM 
is still monotonically decreasing. The other parameters are: 
A = 5 cm” 1 , 7 = 50 cm” 1 , and T = 200 K. (b) The cor¬ 
responding non-Markovianity versus temperature. The re¬ 
sult of TC2 shows finite non-Markovianity, while that from 
TL2 shows very small non-Markovianity and only at low tem¬ 
peratures. The result from HEOM is Markovian due to its 
monotonically-decreasing concurrence. 


again showing a better agreement with the HEOM exact 
results. 

In Fig. 6 (a), 7 is further reduced to 20 cm ” 1 ( 7” 1 = 
265 fs) to investigate the effect of slow environments. As 
the spectral density function is narrower, the correlation 
time of the environment becomes long compared with 
the characteristic time of the system dynamics. Hence 
the information on the system dynamics is more likely 
to be retained in the environment and flow back into the 
system. This back-flow of information in turn affects the 
behavior of the system and results in beating in the con¬ 
currence curves for all methods. As shown in Fig. 6 (b), 
the TC2 model predicts a non-Markovianity much larger 
than the exact results. On the other hand, the TL2 model 
predicts a non-Markovianity in excellent agreement with 
the HEOM results, with only a small under-estimation 














FIG. 6 . (a) The concurrence versus time for TC2 (dot- 

dashed), HEOM (dashed), and TL2 (solid). The 7 value is 
further reduced to 20 cm ^ 1 ( 7 _1 = 265 fs). The other pa¬ 
rameters are the same as those in Fig. 5(a). The information 
on the system dynamics can possibly flow back from the en¬ 
vironment into the system and in turn leads to wavy concur¬ 
rence curves, (b) The corresponding non-Markovianity values 
versus temperature. These non-Markovianity values increase 
prominently as a result of the reduced 7 value. TC2 shows 
larger non-Markovianity values, while TL2 shows good agree¬ 
ment with the HEOM. 


of the non-Markovianity in this set of parameters. 

The above comparisons exhibit an interesting tendency 
for AfM to increase with A and temperature. Several rel¬ 
evant theoretical and experimental works have reported 
[ 3-2 ] that strong system-environment correlations are 
helpful for maintaining quantum coherence even at high 
temperatures. As a result, higher temperature may in 
turn activate more phonon modes in the environment 
without destroying the quantum coherence significantly. 
This provides more channels via which the system can 
interact with the environment. In the language of quan¬ 
tum information science, smaller 7 and strong system- 
environment correlation may help to preserve the dynam¬ 
ical information; while larger A and higher temperature 
may increase the possibility that this information can 
flow back into the system from the environment. Con¬ 


sequently, this increase of Af Ad with larger temperature 
and A is a result of the competition between the back-flow 
of information and thermal fluctuations. Meanwhile, the 
magnitude of the concurrence is reduced by the stronger 
random fluctuations in the environment. 


V. CONCLUSIONS 

In summary, we first investigate the dynamics of two 
perturbative second-order QME methods, TC2 and TL2, 
and compare their results with the numerically-exact re¬ 
sults calculated by HEOM. We find that TC2 can approx¬ 
imate the HEOM population better than TL2. However, 
a drawback of the TC2 model is its over-estimation of 
the coherence. This drawback results in the TC2 model 
predicting too much beating behavior in the population 
dynamics and limits the accuracy of TC2. In constrast, 
the TL2 model predicts sluggish dynamics and loss of co¬ 
herence faster than that of the exact HEOM. As a result, 
the population tends to reach thermal equilibrium too 
rapidly. 

To further investigate the dynamics and establish a 
benchmark for the performance of perturbative QMEs, 
we combine the concept of Choi-Jamiolkowski isomor¬ 
phism [23, 24], entanglement with an ancilla [43], and a 
measure of non-Markovianity [25] to provide a quantita¬ 
tive way to determine how much the coherence dynamics 
and memory effects are deviating from the exact result. 
This provides a deep physical insight on the effects of 
each parameter and a single quantity to determine how 
much the QME dynamics deviates from the exact results. 
Here we find that the non-Markovian measure indicates 
that the TL2 approximates HEOM better than TC2 in 
terms of the coherence dynamics and memory effects for 
the dimer system studied here. In addition, while it is 
well understood that the reorganization energy A and 
temperature enhance the effect of thermal fluctuations 
in the environment on the system, increasing these pa¬ 
rameters can have surprising results. In particular, our 
results show that higher temperature increases informa¬ 
tion back-flow from the environment, thus increasing the 
non-Markovianity of the system dynamics, even though 
the concurrence itself undergoes faster decay. These re¬ 
sults could have important implications in the theoretical 
modeling of electronic coherence in photosynthetic sys¬ 
tems [1, 2, 31]. 
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APPENDIX: FULL EXPRESSIONS FOR THE TC2 AND TL2 QUANTUM MASTER EQUATIONS 

The detailed expression of the TC2 integrodifferential QME Eq. (10) is given by 
d 

_ [ e _ ““'' 3 V«,/3(f)] = -iu<x,pexp[-iu a ,pt] p a ,p(t) 


+ / G(t - p^(r)dT 

w J ° 

- V A a , v A v ^ f G(t - r ) e ^-(‘-O e --„^^ i/3 ( T ) dT 

+ V A atV A„ iP f G* (t - ry^^-^e-^p^dT 

w Jo 

- V A^A^p f G*(t - p a y T )dT, 

J 0 


(17) 


where G(t) is the correlation function defined by Eq. (6). Whereas the detailed expression of TL2 QME in Eq. (12) 
is given by 

d 

_ [e~ Vj}a ' llt p aiP {t)\ = -iMa,pexp[-iu at pt]p ai p(t) 



[ G(t - ry^^-^e 

Jo 

/. L,V 

/*£ 

— V A A 

/ G(t - r)e^- (t - T) e- 

/o 


rt 


/ G*{t-T)e iu ^ t ~ T \ 
Jo 

fi,iy 

ft 


/ G^t-ry^eAt-Q, 

[1,1/ 

/ 0 


e Wa ' , ' t p a yt)dr. 


(18) 
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